"""
1.  读取bdf文件, 根据eid规则区分精细网格与自然网格    结构化解析
2.  读取op2文件, 获取每个单元的bottom与top的stress结果数组,需要考虑遍历多个工况    结构化解析
3.  根据公式计算    分布式算法
    3.1 应力坐标转换（正应力 切应力）
    3.2 遍历10°~170°
    3.3 获取36对应力数值(上下表面 * 每个角度)
    3.4 某个单元在所有工况下的   某个表面某个方向的应力最大与最小值  eid_num * 2 * 18 个数组
4. 外部数据,依据mid   数据集成引擎
    4.1 材料的CFQ
    4.2 材料的q
    4.3 材料的S,插值算法
5. 根据疲劳循环次数N判断合格
"""

from pyNastran.op2.op2 import read_op2
from pyNastran.op2.op2 import OP2
from pyNastran.op2.tables.oes_stressStrain.real.oes_plates import RealPlateStressArray
from pyNastran.bdf.bdf import BDF
from pyNastran.bdf.cards.elements.shell import CTRIA3, CQUAD4
import numpy as np
import re
from re import search, I
import pandas as pd
from scipy.interpolate import interp1d


# CFQcutoff数据
cfq_path = r"D:\SF-CAE\疲劳\疲劳\FatigueTaskCombine\config\fatigue\CFQcutoff.xlsx"
cfq_df = pd.read_excel(cfq_path)
cfq_dict = {}
for index, row in cfq_df.iterrows():
    cfq_dict[row.iloc[4]] = row.iloc[-1]


def get_mat_prop(mid: int) -> dict:
    """根据mid获取相应的材料性能数据

    Args:
        mid (int): 材料编号

    Returns:
        dict: 材料性能数据
    """
    if mid not in cfq_dict:
        # print(f"{mid} 不存在材料数据库中,默认不存在疲劳风险,忽略")
        return
    else:
        prop = {}
        if re.match(r"^3[0-9]{5}$", str(mid)):
            # print("这是钛")
            prop["name"] = "ti"
            prop["C"] = 0.9
            prop["fm0"] = 620
            prop["q"] = 6.5
        if re.match(r"^4[0-9]{5}$", str(mid)):
            # print("这是钢")
            prop["name"] = "st"
            prop["C"] = 0.91
            prop["fm0"] = 930
            prop["q"] = 7.5
        else:
            # print("这是铝")
            prop["name"] = "al"
            prop["C"] = 0.86
            prop["fm0"] = 310
            prop["q"] = 4.0
        prop["cutoff"] = cfq_dict[mid]
        return prop


class ShellStressItem:

    fiber_distance = 0
    oxx = 1
    oyy = 2
    txy = 3
    angle = 4
    omax = 5
    omin = 6
    von_mises = 7


def check_eid(eids: list[int]):
    """区分网格eids 自然/精细"""
    refine_reference_eids = [
        list(range(10000001, 15000000)),
        list(range(15000001, 20000000)),
        list(range(20000001, 29999999)),
        list(range(30000001, 36999999)),
        list(range(37000001, 39999999)),
        list(range(40000001, 44999999)),
        list(range(50000001, 54999999)),
        list(range(60000001, 64999999)),
        list(range(70000001, 79999999)),
        list(range(80000001, 89999999)),
        list(range(90000001, 94999999)),
        list(range(91000001, 91499999)),
        list(range(91500001, 91999999)),
        list(range(92000001, 92999999)),
        list(range(93000001, 93699999)),
        list(range(93700001, 93999999)),
        list(range(94000001, 94499999)),
        list(range(95000001, 95499999)),
        list(range(96000001, 96499999)),
        list(range(97000001, 97999999)),
        list(range(98000001, 98999999)),
        list(range(99000001, 99499999)),
    ]
    set_refine_eids = set(item for subset in refine_reference_eids for item in subset)
    rough_eids = [item for item in eids if item not in set_refine_eids]
    refine_eids = [item for item in eids if item in set_refine_eids]
    # 前者代表自然网格id,后者代表精细网格id
    return rough_eids, refine_eids


# 读取bdf模型文件
bdf_path = r"D:\SF-CAE\疲劳\疲劳\疲劳\model\C919_S32_D20200519.bdf"
bdf_reader = BDF()
bdf_reader.read_bdf(bdf_path, punch=True, validate=False)

# 自然网格/精细网格
rough_eids, refine_eids = check_eid(bdf_reader.element_ids)


# 读取op2结果文件
op2_path_list = [
    r"D:\SF-CAE\疲劳\疲劳\疲劳\result\C919_S32_14160402.op2",
    r"D:\SF-CAE\疲劳\疲劳\疲劳\result\C919_S32_14160403.op2",
]
op2_reader_list = [read_op2(path) for path in op2_path_list]
subcase_lists = [list(reader.subcase_key.keys()) for reader in op2_reader_list]
subcase_ids = [item for sub_list in subcase_lists for item in sub_list]

def get_card_by_eid(eid: int, bdf_reader: BDF) -> str:
    """根据eid获取单元类型名,目前只判断ctrai3和cquad4

    Args:
        eid (int): 单元编号id
        bdf_rader (BDF): bdf读取对象

    Raises:
        f: 不支持其他类型

    Returns:
        str: 单元类型字符串
    """
    ele = bdf_reader.Element(eid)
    if isinstance(ele, CTRIA3):
        return "ctria3"
    elif isinstance(ele, CQUAD4):
        return "cquad4"
    else:
        raise f"str({type(ele)}) is not supported"


def get_row_index_by_eid_location(
    stress: RealPlateStressArray, eid: int, location: int
) -> int:
    """根据指定的eid和location(0代表bottom 1代表top)获取到对应的行索引

    Args:
        stress (RealPlateStressArray): _description_
    """
    if location not in (0, 1):
        raise "value of location must be equaled 0 or 1"
    ele_node = stress.element_node
    for index, item in enumerate(ele_node):
        if (not item[1]) and (item[0] == eid):
            if location:
                return index + 1
            else:
                return index


def get_stress_by_eid_location_item(
    eid: int, location: int, item_index: list[int], bdf_reader: BDF, op2_reader: OP2
) -> list:
    """根据指定的单元编号、应力位置、应力项索引获取对应所有工况的结果

    Args:
        eid (int): 单元编号
        location (int): 应力位置(0代表bottom 1代表top)
        item_index (int): 应力项索引
        bdf_reader (BDF): bdf读取对象
        op2_reader (OP2): op2读取对象

    Returns:
        list: 结果列表(个数为工况数)
    """
    for index in item_index:
        if index not in list(range(0, 8)):
            raise "value of item_index must be an integer between 0 and 7(included)"
    # 单元类型
    card = get_card_by_eid(eid, bdf_reader)
    # 工况编号
    subcase_keys = list(op2_reader.subcase_key.keys())
    # 应力
    stress_path = f"stress.{card}_stress"
    stresses = []
    for subcase_key in subcase_keys:
        stress = op2_reader.get_result(stress_path)[subcase_key]
        row_index = get_row_index_by_eid_location(stress, eid, location)
        result = [stress.data[0][row_index][index] for index in item_index]
        stresses.append(result)
    return stresses


def get_axis_stress(stress_list: list[list]) -> list[float]:
    result = []
    for stress in stress_list:
        result.append(stress[0] - 1 / 3 * stress[1])
    return result


def translate_stress_matrix(stress: list[float], angle: np.float64):
    stress_matrix = np.array([[stress[0], stress[2]], [stress[2], stress[1]]])
    transform_angle_matrix = np.array(
        [[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]
    )
    left_result = np.matmul(np.transpose(transform_angle_matrix), stress_matrix)
    right_result = np.matmul(left_result, transform_angle_matrix)
    transformed_stress_x = right_result[0][0]
    transformed_stress_y = right_result[1][1]
    transformed_single_axis_stress = transformed_stress_x - 1 / 3 * transformed_stress_y
    return transformed_single_axis_stress


def modify_ratio(ratio: float):
    return ((1 - ratio) / 0.94) ** 0.58


def get_max_equal_stress(ratio: float, stress_max: float, mat_prop: dict):
    if ratio <= 0.6:
        return modify_ratio(ratio) * stress_max
    else:
        num1 = 2.5 * (1 - ratio) * mat_prop["fm0"] * stress_max
        num2 = (
            mat_prop["fm0"]
            + 2 * (1 - ratio) * stress_max
            - 0.5 * (1 + ratio) * stress_max
        )
        return modify_ratio(0.6) * num1 / num2


# 根据单元编号区分所属机翼还是机身
class RegexMap:
    def __init__(self, n_dic, val):
        self._items = n_dic
        self.__val = val

    def __getitem__(self, key):
        key = str(key)
        for regex in self._items.keys():
            if search(regex, key, I):
                return self._items[regex]
        return self.__val


def get_section_name_from_eid(eid: int):
    # 精细模型单元编号，共8位
    section_dict = {
        r"^1[0-4][0-9]{6}$": "S11",
        r"^1[5-9][0-9]{6}$": "S12",
        r"^2[0-9][0-9]{6}$": "S20",
        r"^3[0-6][0-9]{6}$": "S31",
        r"^3[7-9][0-9]{6}$": "S32",
        r"^4[0-4][0-9]{6}$": "S40",
        r"^5[0-4][0-9]{6}$": "S50",
        r"^6[0-4][0-9]{6}$": "S60",
        r"^7[0-9][0-9]{6}$": "S70",
        r"^8[0-9][0-9]{6}$": "S80",
        r"^9[0-4][0-9]{6}$": "S90",
    }
    reg_section_dict = RegexMap(section_dict, "S00")
    return reg_section_dict[eid]


oef_x_path = r"D:\SF-CAE\疲劳\疲劳\MSCoef_x.xlsx"
oef_z_path = r"D:\SF-CAE\疲劳\疲劳\MSCoef_z.xlsx"
df_oef_x = pd.read_excel(oef_x_path)
df_oef_z = pd.read_excel(oef_z_path)


def get_fatigue_damage_ratio(eid: int, center: list[float]):
    section_num = int(get_section_name_from_eid(eid)[1:])
    if section_num > 40:
        # z向插值
        f = interp1d(df_oef_z.iloc[:, 0], df_oef_z.iloc[:, 1])
        return float(f(center[2])[0])
    else:
        # x向插值
        f = interp1d(df_oef_x.iloc[:, 0], df_oef_x.iloc[:, 1])
        return float(f(center[0]))


# 计算疲劳循环次数
def get_fatigue_cycle(stress_range: list[float], mid: int, S: float):
    # CFQ和q值,如果未找到,忽略
    prop = get_mat_prop(mid)
    if prop:
        cfq = prop["cutoff"]
        q = prop["q"]
        # 根据应力幅计算应力比 应力比对疲劳计算是有影响的
        stress_ratio = min(stress_range) / max(stress_range)
        if stress_ratio >= -2:
            cw = 1
        elif stress_ratio < -10:
            cw = 0.5
        else:
            cw = 0.00893 * stress_ratio**2 + 0.16966 * stress_ratio + 1.3036
        # 应力比修正
        stress_ratio = stress_ratio * cw
        # 修正系数,取保守值1
        f_sd = 1
        max_equal_stress = get_max_equal_stress(stress_ratio, max(stress_range), prop)
        # 根据公式计算疲劳循环次数
        cycles = S * 100000 * (cfq / max_equal_stress) ** q * f_sd / 10**10
        return cycles
    else:
        # 返回一个极大值
        # print(f"{mid} 没有被找到!!!!!!!!!!")
        return 99999999


print("流程开始")

# for eid in refine_eids:

from concurrent.futures import ThreadPoolExecutor, as_completed
import threading

import queue

# danger_eids = queue.Queue()
# lock = threading.Lock()
completed_tasks = 0
total_tasks = int(100000 * 4.1)
# progress_threshold = 0.1


import sys
import os
sys.path.append(
    os.path.dirname(os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
)
from bdf_transformer_2.bdf_transformer import BDFTransformer

def get_assembly(bdf_dir):
    index = 0
    assembly_info = []
    # comp与cid的映射关系
    map_comp_info = {}
    # pname与pid的映射
    map_prop_info = {}
    for root, _, files in os.walk(bdf_dir):
        for file in files:
            if not file.endswith(".bdf"):
                pass
            else:
                bdf_transformer = BDFTransformer(os.path.join(root, file))
                for info in bdf_transformer.get_comps_info():
                    map_comp_info[info["cid"]] = info["cname"]
                    map_prop_info[info["pid"]] = info["pname"]
                if bdf_transformer.file_category == "mesh":
                    index += 1
                    assembly_info.append(
                        {os.path.splitext(file)[0]: bdf_transformer.get_comp_eids_ref()}
                    )
    return assembly_info, map_comp_info, map_prop_info

assembly = get_assembly(os.path.dirname(bdf_path))[0]
summary_info = []

header = ["G70 Number", "ELEMID", "LF", "Params", "CFQ", "fm0", "s"]
header.extend(subcase_ids)


def evaluate_fatigue(eid):
    result = {"bottom": {}, "top": {}}
    ele = bdf_reader.Element(eid)
    center = ele.Centroid()
    # 材料id
    mid = ele.material_ids[0]

    if get_mat_prop(mid) is not None:
        # 单元类型
        card = get_card_by_eid(eid, bdf_reader)
        # 单元的默认σx σy τxy
        # 单个eid在 所有工况 bottom/up面 的应力结果
        stress_bottom = []
        stress_top = []
        for op2_reader in op2_reader_list:
            # 下表面
            stress_bottom.extend(
                get_stress_by_eid_location_item(
                    eid,
                    0,
                    [ShellStressItem.oxx, ShellStressItem.oyy, ShellStressItem.txy],
                    bdf_reader,
                    op2_reader,
                )
            )
            # 上表面
            stress_top.extend(
                get_stress_by_eid_location_item(
                    eid,
                    1,
                    [ShellStressItem.oxx, ShellStressItem.oyy, ShellStressItem.txy],
                    bdf_reader,
                    op2_reader,
                )
            )
        S = get_fatigue_damage_ratio(eid, center)
        # 处理所有角度
        axis_stress_bottom = []
        axis_stress_top = []
        for i in range(0, 18):
            angle = np.pi / 18 * i
            temp_stress1, temp_stress2 = [], []
            for stress in stress_bottom:
                temp_stress1.append(translate_stress_matrix(stress, angle))
            for stress in stress_top:
                temp_stress2.append(translate_stress_matrix(stress, angle))
            # 第一个元素为最小单轴应力  第二个元素为最大单轴应力  本质上是计算应力幅值
            axis_stress_bottom.append(temp_stress1)
            axis_stress_top.append(temp_stress2)
            result["bottom"][i] = get_fatigue_cycle(
                [min(temp_stress1), max(temp_stress1)], mid, S
            )
            result["top"][i] = get_fatigue_cycle(
                [min(temp_stress2), max(temp_stress2)], mid, S
            )


        # 判断该单元的疲劳循环次数是否满足要求
        all_cycles = list(result["bottom"].values()) + list(result["top"].values())
        danger_angle_index_bottom, danger_angle_index_top = [], []
        for index, cycle in enumerate(list(result["bottom"].values())):
            if cycle <= 10**5:
                danger_angle_index_bottom.append(index)
        for index, cycle in enumerate(list(result["top"].values())):
            if cycle <= 10**5:
                danger_angle_index_top.append(index)
        if len(danger_angle_index_top) + len(danger_angle_index_bottom) == 0:
            pass
        else:
            for info in assembly:
                for k, v in info.items():
                    for cid, eids in v.items():
                        if eid in eids:
                            g70 = k + ".bdf"
                            break
                    break
            CFQ = get_mat_prop(mid)["cutoff"]
            fm0 = get_mat_prop(mid)["fm0"]
            s = get_mat_prop(mid)["q"]
            for index in danger_angle_index_bottom:
                if index < 10:
                    elemid = f"{eid}@10{index}"
                else:
                    elemid = f"{eid}@1{index}"
                base = [g70, elemid, "", "1.0;1.5", CFQ, fm0, s]
                base.extend(axis_stress_bottom[index])
                summary_info.append(base)
            for index in danger_angle_index_top:
                if index < 10:
                    elemid = f"{eid}@20{index}"
                else:
                    elemid = f"{eid}@2{index}"
                base = [g70, elemid, "", "1.0;1.5", CFQ, fm0, s]
                base.extend(axis_stress_top[index])
                summary_info.append(base)



import time

start = time.time()


progress_points = [i * total_tasks // 10 for i in range(1, 11)]

with ThreadPoolExecutor(max_workers=16) as executor:
    futures = {executor.submit(evaluate_fatigue, eid): eid for eid in refine_eids[:total_tasks]}

    completed_tasks = 0
    last_printed_progress = -1

    # 等待任务完成
    for future in as_completed(futures):
        completed_tasks += 1
        current_progress = (completed_tasks * 100) // total_tasks

        # 检查是否达到新的10%进度点
        if current_progress // 10 > last_printed_progress // 10:
            print(f"Progress: {current_progress}%")
            last_printed_progress = current_progress

end = time.time()

print(f"计算耗时 {end-start}")

df = pd.DataFrame(summary_info, columns=header)
df.to_csv("out.csv", index=False)